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The  limiting  current  density  is  an  important  characteristic  quantity  in  solid  oxide  fuel  cells  (SOFCs).  High 
concentration  overpotential  is  often  used  to  explain  the  limiting  current  density  assuming  a  high  tortu¬ 
osity  or  limited  surface  diffusion  in  the  vicinity  of  the  three-phase  boundary.  Most  membrane-electrode 
assembly  models  of  SOFC  fail  to  predict  the  limiting  current  density,  even  for  hydrogen,  when  using  physi¬ 
cally  reasonable  values  of  tortuosity  and  considering  the  short  residence  time  of  the  adsorbed  species  near 
the  three-phase  boundary.  In  this  paper,  a  one-dimensional  model  for  the  transport-chemistry  interac¬ 
tions  in  SOFCs  is  described.  The  model  is  based  on  a  comprehensive  approach  that  includes  the  dusty-gas 
model  for  gas  transport  in  the  porous  electrodes,  detailed  heterogeneous  elementary  reaction  kinetics 
for  the  thermo-chemistry  in  the  anode,  and  detailed  electrode  kinetics  for  the  electrochemistry.  Correct 
values  for  the  Knudsen  diffusion  coefficients  are  used.  We  apply  the  unsteady  form  of  the  conservation 
equations,  allowing  for  the  analysis  of  the  response  of  the  cell  to  external  dynamics.  Results  of  our  model 
are  compared  with  experimental  data,  showing  good  agreement  over  a  wide  range  of  the  current  density, 
but  fail  to  predict  the  limiting  current  density  accurately  when  the  hydroxyl  oxidation  charge-transfer 
reaction  is  assumed  to  be  the  rate  limiting  reaction.  To  obtain  accurate  predictions  of  the  limiting  current 
density,  we  analyze  the  possibility  that  different  steps  can  be  rate-limiting  reactions  in  the  electrochem¬ 
istry  model  of  hydrogen  oxidation.  We  use  recent  measurements  on  the  three-phase  boundary  area  and 
take  into  account  the  surface  diffusion  and  competitive  adsorption  to  determine  possible  rate  limiting 
reactions  at  high  current  density.  We  show  that  a  rate-limiting  switchover  model,  in  which  the  reaction 
limiting  the  overall  kinetic  rate  becomes  the  hydrogen  adsorption  at  the  anode,  may  be  required  to  explain 
the  experimentally  measured  limiting  current  density  over  a  range  of  operating  conditions. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  are  considered  attractive  alternatives  to  combustion 
engines  because  of  their  high  theoretical  efficiency.  Solid  oxide  fuel 
cells  (SOFCs)  have  special  advantages  over  low  temperature  fuel 
cells  such  as  polymer  electrolyte  membrane  fuel  cells  because  they 
exhibit  high  conversion  efficiency  at  high  temperatures  while  con¬ 
suming  hydrocarbons  directly.  Current  SOFC  research  focus  areas 
include  reducing  the  operating  temperature  and  eliminating  the 
need  for  reformers  while  utilizing  hydrocarbon  fuels.  In  order  to 
achieve  these  objectives  while  improving  the  finite  current  effi¬ 
ciency  and  optimizing  the  design  of  SOFCs,  mathematical  models 
are  indispensable. 

The  limiting  current  density  is  an  important  characteristic  quan¬ 
tity  in  SOFCs.  Efforts  to  predict  the  limiting  current  density  have  not 
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been  successful.  The  high  concentration  overpotential,  determined 
by  the  reactants  and  products  transport  limitations  that  reduce 
the  reactants’  concentrations  at  the  three-phase  boundary  (TPB), 
is  often  used  to  explain  the  limiting  current  density.  The  concen¬ 
tration  overpotential  depends  on  the  transport  properties  of  the 
electrodes  including  their  porosities  and  tortuosities.  To  predict 
the  measured  limiting  current  density,  some  membrane-electrode 
assembly  (MEA)  models  assume  anode  tortuosities  in  the  range 
of  10-17.  These  values  are  too  high,  considering  that  the  reported 
range  for  porous  sintered  ceramics  is  2-10,  and  most  often  below  6 
[  1  ].  To  avoid  making  this  assumption,  Williford  et  al.  [  1  ]  introduced 
surface  diffusion  effects  into  Fields  diffusion  model,  by  adjusting 
the  diffusion  coefficients  assuming  that  the  competitive  adsorption 
and  surface  diffusion  contribute  to  the  concentration  overpotential. 
The  adjusted  diffusion  coefficient  was  introduced  by  combining  a 
gas  diffusion  coefficient  and  a  surface  diffusion  coefficient,  and  was 
applied  to  the  gas  species  transport  throughout  the  entire  thickness 
of  electrodes.  This  approach  overestimated  the  diffusion  resistance 
because:  (i)  the  diffusion  coefficient  adjusted  using  the  surface 
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diffusion  coefficient  is  applied  to  the  gas  species  transport;  and 
(ii)  the  adjusted  coefficient  is  applied  to  the  whole  electrode  even 
though  surface  diffusion  plays  a  role  only  in  the  vicinity  of  the  TPB, 
within  tens  of  nanometers. 

In  this  paper,  we  propose  other  mechanisms  that  can  explain  the 
cell  behavior  around  the  limiting  current  density.  For  this  purpose, 
we  develop  an  accurate  transport-chemistry  interaction  model  for 
the  MEA,  using  the  dusty-gas  model  (DGM)  to  describe  transport; 
a  detailed  mechanism  for  the  thermochemistry  in  the  anode;  and 
multistep  reaction  mechanisms  to  model  the  charge-transfer  reac¬ 
tions  at  the  anode  and  the  cathode  surfaces.  We  show  that  the 
model  can  predict  the  polarization  curve  accurately  over  most  of  the 
range  of  current  density,  while  using  the  assumption  that  hydroxyl 
ion  oxidation  reaction  is  the  limiting  step  in  the  anode  electrochem¬ 
istry  model.  However,  this  assumption  fails  to  predict  the  limiting 
current  density  accurately.  We  explore  other  rate  limiting  reac¬ 
tions  at  high  current  density  and  show  that  hydrogen  adsorption 
at  the  anode  can  become  sufficiently  slow  to  cause  rapid  increase 
of  the  activation  overpotential  at  high  currents.  Using  a  modified 
approach  to  the  selection  of  the  rate  limiting  reaction,  we  obtain 
better  prediction  of  the  limiting  current  density. 

In  the  following  sections,  a  detailed  model  to  calculate  each 
overpotential  and  the  way  to  implement  them  in  the  simulation 
code  will  be  described.  We  will  validate  the  model  against  available 
experimental  results  and  propose  a  new  rate-limiting  switch-over 
model  to  improve  the  prediction  of  the  limiting  current  density. 
The  paper  ends  with  conclusions  and  suggestions  for  future  work. 


2.  Model  description 


The  objective  of  the  model  is  to  calculate  the  polarization  curve 
of  a  SOFC  over  the  entire  operating  range  of  voltage-current  density, 
and  for  different  fuel  concentrations.  We  adopt  a  one-dimensional 
approach  to  model  transport-chemistry  interactions  within  the 
MEA  and  limit  our  validation  to  button  cell  results.  Extension  to 
multi-dimensional  models  will  be  attempted  in  the  future.  We 
assumed  that  the  temperature  is  constant  and  uniform  through  the 
MEA. 

The  equilibrium  potential,  Erev,  is  calculated  from  the  Gibb’s  free 
energy  of  the  reaction: 
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rev  — 
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where  vt  is  the  stoichiometric  coefficient  of  a  constituent,  3T  is  the 
universal  gas  constant  (J  mol-1  K-1 ),  T  is  the  temperature  (K),  p*  is 
the  partial  pressure  of  gas  species  i  (atm),  and  the  Gibbs  free  energy 
of  reaction,  Ag°  (J  mol-1 )  is  evaluated  for  all  species  at  atmospheric 
pressure.  Erev  is  also  the  open  circuit  voltage  (OCV)  with  no  fuel 
leakage. 

At  finite  current,  the  open  circuit  voltage  is  reduced  by  losses, 
including  the  ohmic  overpotential  associated  with  ion  transport 
through  the  electrolyte  and  electron  transport  through  the  elec¬ 
trodes,  the  activation  overpotentials  associated  with  the  energy 
barriers  of  the  charge-transfer  reactions,  and  the  concentration 
overpotentials  associated  with  the  gas-phase  species  transport 
resistance  through  the  electrodes.  Thus,  the  operating  cell  voltage, 
E ,  can  be  written  as 


Ecell  ~  Erev  ~  P  cone,  a  ~  If  a, a  ~  If  ohm  ~  If  cone, c  ~  Pa,c 


(2) 


2.1  Concentration  overpotential 

At  open-circuit  conditions,  i.e.  zero  current  flow,  the  species  con¬ 
centrations  at  the  TPB  are  the  same  as  those  of  the  bulk  channel’s 
gases.  Under  finite  current  conditions,  species  concentrations  at 
the  TPB  are  different,  because  of  the  finite  transport  flux  across  the 
electrodes.  While  evaluating  the  actual  electrochemical  potential 
of  the  fuel  cell,  the  relevant  reactants  and  products  concentrations 
are  those  at  the  TPB.  The  potential  difference  associated  with  the 
gas  species  concentration  change  across  the  electrodes  is  the  con¬ 
centration  overpotential: 

RT 

Pconc  =  [£rev]at  the  channel  —  [^revlatthe  TPB  =  ~ 
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Next  we  develop  a  model  for  computing  the  concentrations  of  gas 
species  at  the  TPB. 

2.1.1.  Conservation  equation 

The  conservation  equation  of  gas-phase  species  in  the  reactive 
porous  media  is 

—  AsSSurf,k  +  $gas,k  ~  ^  'Jk  (k  =  1 ,  .  .  .  ,  Kg)  (4) 

where  ck  is  the  concentration  of  gas  species  k  (molm ~3)Jk  is  the 
molar  flux  of  gas  species  k  (mol  m-2  s-1 ),  ssurjk  is  the  surface  pro¬ 
duction  rates  of  the  gas  species  k  by  the  heterogeneous  reactions 
(mol  m-2  s-1 ),  As  is  the  specific  catalyst  area  per  unit  volume  of  the 
electrode  (m-1 ),  sgask  is  the  production  rates  of  the  gas  species  k  by 
the  homogeneous  reactions  (mol  m-3  s-1 ),  and  I<g  is  the  total  num¬ 
ber  of  gas  species.  The  molar  flux  can  be  determined  by  the  Fields 
model  (FM)  or  the  DGM.  The  production  rates  of  the  gas  species  are 
obtained  from  the  thermo-chemistry  model. 

The  surface  species  conservation  equation  is  as  follows: 

(k=\,...,Ks)  (5) 

where  csurfk  is  the  concentration  of  surface  species  k  (molm-2), 
sk  is  the  production  rate  of  the  surface  species  k  by  the  hetero¬ 
geneous  reactions  (molm-2s-1),  and  I<s  is  the  total  number  of 
surface  species.  Unlike  gaseous  species,  surface  species  are  effec¬ 
tively  immobile  on  length  scales  larger  than  an  individual  catalyst 
particle.  Hence,  the  surface  species  transport  over  macroscopic  dis¬ 
tance  is  assumed  negligible  [2]. 


2.1.2.  Transport 

The  fluxes  Jk  are  computed  using  the  DGM  [3].  DGM  incorpo¬ 
rate  physical  phenomena  beyond  those  described  by  FM,  such  as 
osmotic  diffusion,  reverse  diffusion,  and  diffusion  barrier  [4]: 


B°  -^-Vp 

//„«  d?m  p 


(6) 


where  ct  =  p/9TT  is  total  molar  concentration  (molm-3). 

The  effective  Knudsen  diffusion  coefficient  for  the  component  i, 
DeM,  in  the  multicomopnent  mixture  gas  is  governed  by  [3]: 


where  pConc,a  and  pC0nCiC  are  the  concentration  overpotentials  at  the 
anode  and  the  cathode,  pa,a  and  r)a,c  are  the  corresponding  acti¬ 
vation  overpotentials  at  the  anode  and  the  cathode,  and  pohm  is 
the  total  ohmic  overpotential.  Models  for  these  overpotentials  are 
developed  next. 


e  _  do  £  j  8RT 
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where  M*  is  the  molar  mass  of  species  i  (kg  mol-1 ),  D|,  is  the  effec¬ 
tive  binary  diffusion  coefficient  in  the  porous  medium,  and  D?.  is 
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related  to  the  corresponding  ordinary  binary  diffusion  coefficient 
Dy  as  [4]: 


where  e  is  the  porosity  and  r  is  the  tortuosity.  The  porosity  is 
defined  as 

s=^mL_  (9) 

v  material 

in  which  Vvold  is  a  void  volume  and  V material  is  the  superficial  volume 
of  the  material. 

The  tortuosity  is  defined  as 


r  = 


(10) 


where  le  is  the  effective  length  between  two  points  through  the 
pores  and  l  is  the  straight  distance  between  the  same  two  points. 

According  to  Champan-Enskog  kinetic  theory,  the  binary  diffu¬ 
sion  coefficient  Dy  is  [5]: 


D„  =  5.8765  x  10“9 
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where  T2Dy  is  the  dimensionless  collision  integral  function  [5]: 
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2.1.3.  Thermo -chemistry 

Because  of  the  relatively  high  operating  temperatures  and 
the  catalytic  surfaces  in  the  anode  structure,  various  ther¬ 
mochemical  reactions  occur  within  the  anode,  such  as  steam 
reforming,  water-gas  shift,  partial  oxidation,  and  carbon  forma¬ 
tion.  Thermochemistry  has  usually  been  handled  using  significantly 
simplifying  assumptions,  such  as  local  equilibration  of  reforming 
and  water-gas-shift  chemistry  [1,7],  or  global  reaction  kinetics  [8]. 
Recently,  detailed  heterogeneous  elementary  kinetics  models  of 
methane  on  a  Ni  surface  have  been  established  and  validated  over 
a  wide  range  of  conditions  [9]. 

Gas-phase  chemistry  within  the  pores  is  neglected  because  the 
heterogeneous  thermochemical  reactions  are  considerably  faster 
than  the  homogeneous  thermochemistry  and  the  probability  for 
gas-gas  collisions  is  low  when  the  pore  space  is  comparable  to  the 
mean  free-path  length  [9].  Thus, 

W  =  °  ik  =  \,...,Kg)  (20) 

The  surface  mechanism  of  methane  reforming  and  oxidation  over 
nickel  has  been  suggested  [10].  The  mechanism  was  initially  devel¬ 
oped  and  validated  using  Ni-coated  honeycomb  monoliths  for 
the  temperature  range  from  700  to  1300  K.  The  reaction  mecha¬ 
nism  consists  of  6  pairs  of  adsorption/desorption  reactions  for  6 
gas  species  and  15  pairs  of  surface  reactions  among  12  adsorbed 
species.  We  included  all  44  heterogeneous  elementary  reactions 
in  the  model.  Details  on  how  to  apply  detailed  thermochemical 
reforming  to  our  model  can  be  found  in  [11  ]. 


where  kB  is  the  Boltzmann  constant  (J  K-1 )  and  Sy  is  the  character¬ 
istic  Lennard-Jones  energy  (J).  Here,  cry  and  Sy  are  calculated  from 
the  individual  parameters  using  the  approximate  equations  [5]: 


The  mixture  viscosity,  iivmix,  is  [5]: 
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The  viscosity  of  each  species  is  determined  by  [5]: 


Hv  =  8.4411  x  10“ 
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2.2.  Activation  overpotential 

Charge-transfer  processes  are  among  the  least  understood 
aspects  of  fuel-cell  chemistry.  Analytical  approaches  using  a  single 
global  charge-transfer  reaction  have  often  been  used  to  describe 
electrochemical  kinetics,  leading  to  the  Butler-Volmer  equation 
[1,12,13].  There  have  been  some  efforts  to  develop  detailed  charge- 
transfer  kinetics  in  the  form  of  elementary  reaction  steps,  in 
a  manner  that  resembles  the  treatment  of  heterogeneous  ther¬ 
mochemical  kinetics.  However,  a  clear  understanding  of  the 
electrode  kinetics  does  not  exist  yet.  Regarding  the  anode,  there 
are  several  outstanding  issues.  According  to  the  literature,  adsorp¬ 
tion/desorption,  surface  diffusion,  the  formation  of  hydroxyl  and 
a  charge-transfer  reaction  are  feasible  rate-limiting  reaction  steps 
in  a  simplified  SOFC  anode  model.  It  is  not  evident  whether  the 
chemical  and  the  electrochemical  reactions  occur  only  on  the  sur¬ 
faces  of  Nj  and  of  yttria-stabilized  zirconia  (YSZ),  or  whether  the 
bulk  material  is  also  active. 

Zhu  et  al.  [9]  applied  a  detailed  electrode  kinetics  model  and 
obtained  a  Butler-Volmer  formalism  using  a  rate-limiting  step.  This 
approach  provides  qualitative  information  about  some  important 
functional  dependencies  such  as  the  reaction  order  in  the  exchange 
current  density,  and  enables  comparison  with  experimental  results 
that  have  been  interpreted  using  parameters  in  the  Butler-Volmer 
equation.  In  the  following,  we  follow  the  same  approach. 


The  permeability  B0  is  an  experimentally  determined  characteristic 
parameter  of  the  porous  matrix  structure  including  the  porosity  and 
tortuosity  factors.  If  an  electrode  is  assumed  to  be  an  aggregated 
bed  of  spherical  particles  with  diameter  dp  (m),  the  permeability 
can  be  expressed  by  the  Kozeny-Carman  relationship  [6]: 


Bo 


dj  s 3 

180  (1  _  £)2 


(19) 


The  transport  of  gaseous  species  through  porous  electrodes  is 
affected  by  the  microstructure  of  the  electrodes,  particularly,  the 
porosity,  permeability,  pore  size,  and  tortuosity  factor. 


2.2.1.  Anode  activation  overpotential 

We  use  the  five  elementary  reaction  mechanisms  proposed  by 
de  Boer  [14].  In  this  model,  hydrogen  is  adsorbed  only  on  the  nickel 
surface  (Ni)  and  other  surface  species  reside  on  the  electrolyte 
surface  (YSZ).  These  reactions  include  interactions  among  (i)  an 
adsorbed  atomic  hydrogen,  H(Ni),  an  empty  surface  site,  (Ni),  and 
an  electron,  e-(Ni),  on  the  Ni  anode  surface;  (ii)  a  lattice  oxygen, 
Oq  (YSZ),  and  an  oxygen  vacancy,  V0##  (YSZ)  within  the  YSZ  elec¬ 
trolyte;  and  (iii)  hydroxyl  ion  OH_(YSZ),  water  H20(YSZ),  oxygen 
ion  02_(YSZ),  and  empty  sites  (YSZ)  on  the  YSZ  surface.  The  reaction 
mechanism  is 
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Table  1 

The  different  terms  and  parameters  in  BV  equation  when  assuming  that  Reaction  (Rl),  (R2),  (R3)  or  (R4)  is  rate-limiting. 


Adsorption  on  the  Ni  surface: 

H2(g)  +  2(Ni)  ^  2H(Ni)  (Rl) 

Charge-transfer  reactions  at  the  TPB  region: 


other  reaction  rates  match  the  expected  values  of  the  current  den¬ 
sity.  In  the  following,  we  calculate  the  hydrogen  adsorption  rate 
based  on  the  estimate  of  the  TPB  area  using  recent  3D  reconstruc¬ 
tion  of  the  anode/electrolyte  interface. 


H(Ni)  +  02-(YSZ)  4^  (Ni)  +  OH-(YSZ)  +  e”(Ni)  (R2) 

H(Ni)  +  OH-(YSZ)  (Ni)  +  H20(YSZ)  +  e”(Ni)  (R3) 

•  Adsorption/desorption  on  the  YSZ  surface: 

H20(YSZ)  4*  H20(g)  +  (YSZ)  (R4) 

•  Transfer  of  oxygen  ion  between  the  surface  and  the  bulk  YSZ: 

Oq(YSZ)  +  (YSZ)  4*  02-(YSZ)  +  Vq##(YSZ)  (R5) 


The  calculated  residence  time  of  H  on  Ni  surface  is  8  ns  and 
the  corresponding  diffusion  length  is  19  nm  at  750  °C.  Because  the 
adsorbed  hydrogen  desorbs  after  this  small  diffusion  length,  the 
surface  coverage  is  mainly  determined  by  adsorption/desorption. 
Even  within  the  diffusion  length  scale,  50  nm,  away  from  the  TPB 
the  surface  coverage  changes  by  only  0.1%.  Thus,  it  is  reasonable  to 
exclude  the  surface  diffusion  of  hydrogen  from  this  mechanism. 

The  surface  species  interact  through  five  reaction  equations, 
Reactions  (Rl  )-(R5)  and  the  two  conservation  equations  for  Ni  and 
YSZ  surface.  Invoking  the  assumption  that  one  reaction  step  is  rate- 
limiting,  a  simple  analytical  expression  for  the  i-r\a  relationship  can 
be  derived  in  the  Butler-Volmer  form  (BV)  using  the  anodic  and  the 
cathodic  reaction  coefficient,  ki  a  and  ki  c  of  the  rate-limiting  reac¬ 
tion  and  the  equilibrium  constants  I<j  of  other  reactions.  The  form 
of  the  BV  equation  is 


i  =  i0 


(21) 


where  /3a  and  f3c  are  the  anodic  and  cathodic  charge-transfer  coef¬ 
ficient,  respectively.  The  results  of  these  derivations,  assuming  that 
Reaction  (Rl),  (R2),  (R3)  or  (R4)  may  be  the  rate  limiting  reaction, 
are  shown  in  Table  1.  Note  that  different  expressions  are  obtained 
for  the  exchange  current  density,  their  dependence  on  the  par¬ 
tial  pressure  of  the  reactants  and  product,  and  the  charge-transfer 
coefficients.  The  form  of  these  expressions  are  important  in  deter¬ 
mining  the  actual  rate  limiting  reaction  under  different  operating 
conditions,  as  shown  next. 

Zhu  et  al.  [9]  derived  the  BV  type  i-r\a  relationship  assuming 
that  Reaction  (R3)  is  rate-limiting,  and  used  that  expression  in  their 
modeling  analysis.  They  based  their  choice  on  the  argument  that 
the  hydrogen  adsorption  rate  is  several  orders  of  magnitude  higher 
than  the  current  density  of  interest  [15,16]  and  hence  cannot  be 
rate-limiting.  In  the  following  we  show  that  this  is  indeed  a  good 
choice  for  the  rate  limiting  reaction  over  a  wide  range  of  current 
density,  but  not  necessarily  at  high  current  density. 

We  analyzed  the  possibilities  that  reactions  other  than  Reaction 
(R3)  can  be  rate-limiting,  especially  at  high  current  density,  and 
examined  the  possibility  that  hydrogen  adsorption  Reaction  (Rl) 
may  be  the  rate-limiting.  For  this  purpose,  we  need  to  show  that 


(a)  TPB  area: 

Recently,  a  3D  reconstruction  of  the  anode  was  obtained  by 
stacking  the  2D  scanning  electron  microscopy  (SEM)  images 
in  3D  space.  From  the  3D  reconstruction,  the  volume-specific 
TBP  length  was  directly  determined.  It  was  found  to  be 
4.28  x  1012  m-2.  For  the  electrochemical  reactions  to  take  place, 
the  TPB  must  be  connected  to  the  rest  of  the  structure.  That 
is,  the  pores  must  be  connected  through  the  surrounding  pore 
network  to  the  fuel  stream,  the  Ni  phase  connected  to  the 
surrounding  Ni  phase  then  ultimately  to  the  current  collector, 
and  the  YSZ  phase  connected  to  the  bulk  YSZ  electrolyte.  It  is 
reported  that  63%  of  TPBs  were  interconnected  [17].  Moreover, 
it  was  argued  that  the  TPB  width  is  in  the  range  of  10-9-10-1°  m 
in  [18]  and  5  x  10-1°  m  in  [1  ].  The  active  depth  of  the  anode  is 
0(10  |xm)  [19].  Based  on  these  values,  the  TPB  area  per  electrode 
area  can  be  calculated  as  follows: 

Atpb  =  0.63  x  (4.28  x  1012)  x  (10  x  10“6)  x  (5  x  10“10) 

A electrode 

=  0.013  (22) 

(b)  Hydrogen  adsorption  rate  at  the  TPB: 

Using  the  sticking  coefficient  of  hydrogen  on  nickel,  the 
hydrogen  adsorption  rate  can  be  computed  from  the  following 
expression  [20]: 

1231 

where  0V  is  a  vacancy  coverage. 

The  concentration  of  hydrogen  in  the  gas  channel  is  in  the 
order  of  several  molrn-3.  If  all  the  adsorbed  hydrogen  reacts, 
the  order  of  magnitude  of  the  resulting  current  density  is  as 
follows: 

i  =  —  x  F  x  R]  ~  O  (1  Aon-2)  (24) 

A electrode 

This  result  shows  that  the  hydrogen  adsorption  rate  and  the  cur¬ 
rent  density  of  interest  are  indeed  comparable.  Thus,  at  high  current 
density,  hydrogen  adsorption  at  the  anode  may  be  the  rate  limiting 
reaction.  We  should  mention,  however,  that  the  current  density 
estimated  in  Eq.  (24),  based  on  the  TPB  area  and  the  adsorption 
rate  of  hydrogen  on  this  area,  may  be  too  high.  For  instance,  since 
Atpb  is  the  combined  area  of  nickel  and  YSZ  surfaces,  the  actual 
area  of  Ni  surface  for  hydrogen  adsorption  may  be  smaller  than  the 
value  given  in  Eq.  (22).  Furthermore,  the  present  model  assumes 
that  only  hydrogen  can  be  adsorbed  on  the  nickel  surfaces.  How¬ 
ever,  other  species  such  as  H2,  H20,  OH,  and  O  can  also  be  adsorbed 
on  nickel  surfaces  [1].  Considering  the  actual  available  Ni  sur¬ 
face  area  of  the  TPB  for  hydrogen  adsorption,  and  the  competitive 
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adsorption  of  other  species,  we  conclude  that  the  hydrogen 
adsorption  rate,  under  the  condition  of  high  current,  can  be  the  rate- 
limiting  reaction.  Thus,  the  case  when  Reaction  (R1 )  is  rate-limiting 
at  high  current  density  is  included  in  our  analysis.  As  shown  next, 
other  considerations  must  be  taken  into  account  in  defining  the  rate 
limiting  reaction. 

Determining  the  actual  rate  limiting  reaction  among  Reactions 
(R1)-(R4),  requires  careful  analysis,  which  is  based  on  several 
observations  and/or  matching  the  form  of  the  corresponding 
Butler-Volmer  equation  to  experimental  data.  Experimental  results 
have  shown  a  number  of  trends  that  are  relevant  to  the  anode 
charge-transfer  reaction. 

(1) The  anodic  charge-transfer  coefficient,  /3a,  is  greater  than 
the  value  often  assumed  when  a  single  rate-limiting  charge- 
transfer  reaction  is  used  over  the  entire  range  of  current 
density;  that  is,  it  is  greater  than  0.5. 

(2)  It  has  been  reported  that  small  amounts  of  water  added  to 
the  fuel  gas  accelerates  the  electrochemical  charge-transfer 
reaction,  known  as  the  catalytic  effect  of  water  [14,16,21  ].  Con¬ 
sidering  the  global  reaction  in  the  anode,  the  partial  pressure 
of  hydrogen,  a  reactant,  is  expected  to  have  a  positive  reaction 
order,  while  the  partial  pressure  of  water,  a  product,  is  expected 
to  have  a  negative  reaction  order.  This  is,  however,  contrary  to 
what  is  observed  experimentally;  as  it  is,  water  promotes  the 
electrochemical  reactions. 

If  Reaction  (R3)  or  (R4)  is  assumed  to  be  the  rate-limiting  pro¬ 
cess: 

(i)  the  anodic  charge-transfer  coefficient  in  the  BV  becomes  1.5  or 
2,  respectively,  when  anodic  and  cathodic  charge-transfer  coef¬ 
ficients  of  the  elementary  charge-transfer  steps  are  assumed  to 
be  0.5,  and 

(ii)  the  reaction  order  of  water  in  the  exchange  current  density 
becomes  positive. 

Thus,  either  reaction  could  be  rate  limiting  and  compatible  with 
the  above  observation.  However,  there  is  a  distinction.  If  Reaction 
(R4)  is  assumed  to  be  the  rate-limiting  reaction,  the  exchange  cur¬ 
rent  density  does  not  depend  on  the  hydrogen  concentration  and 
the  cathodic  charge-transfer  coefficient  reduces  to  zero.  According 
to  de  Boer,  the  anodic  charge-transfer  coefficient  in  most  cases  was 
close  to  1.5  for  cermet  anode.  For  the  cathodic  branch,  the  charge- 
transfer  coefficient  is  approximately  0.5-1.0.  Furthermore,  for  Ni 
patterned  anodes,  a  considerable  influence  was  found  for  pH2  on 
the  polarization  resistance,  especially  at  a  very  low  partial  pressure 
of  hydrogen  [21  ].  Thus,  we  conclude  that  the  analytic  expression  of 
the  Butler-Volmer  form  exhibits  a  better  agreement  with  experi¬ 
ments  when  Reaction  (R3)  rather  than  Reaction  (R4)  is  assumed  to 
be  the  rate-liming  reaction. 


•  Charge-transfer  and  incorporation  at  the  TPB: 

Oad(c)  +  V0 -  (el)  +  2 e-(c)  O *(e()  +  (c)  (R7) 

In  these  reactions  Oad(c)  is  the  adsorbed  atomic  oxygen  on  the 
cathode  surface  and  (c)  is  an  unoccupied  cathode  surface  site. 


Assuming  that  the  charge-transfer  coefficients  are  0.5,  the  BV 
form  of  i-relationship  simplifies  to 


V29t  T  J 
[,1/2 1,1/2 


2mT  J  J 


‘o  ; 


2lTPBFk^k‘/J(po2K6) 


1/4 


(</2+Pa 


V22) 


By  writing 


1/2. 1/2 


lo2  —  21tpbF1<7  c  k7'  a 

the  exchange  current  density  is  expressed  as 

,  ,*  (Po2/K6)1/4 

i0  =  ic 


(25) 

(26) 

(27) 

(28) 


°21+(Po2/K6)1/2 

For  an  LSM-YSZ  interface,  /C6  is  represented  in  Arrhenius  form  [22] : 


K 6  =  Ao2  exp 


_ ’2_ 

WT 


(29) 


in  which  Ao2  =  4.9  x  108  atm  and  Eq2  =  200kJ/mol. 

Even  though  there  is  a  discrepancy  in  the  exponent  of  oxygen 
partial  pressure  in  the  numerator,  Eq.  (28)  matches  well  the  exper¬ 
imental  results  given  by  [22]: 


*o  =  *o2 


(P02)1/2 

1+(Po2/K6)1/2 


(30) 


The  parameter  i*2  is  taken  here  as  an  adjustable  parameter  that 
is  varied  to  match  the  experimentally  observed  cell  performance 
reported  by  Jiang  and  Virkar  [7]  under  one  set  of  conditions.  It  is 
then  used  unchanged  to  predict  the  cell  performance  under  other 
conditions. 


2.3.  Ohmic  overpotential 

Sources  of  ohmic  losses  in  a  fuel  cell  are  the  ion  flow  resistance 
across  the  electrolyte  and  the  electronic  flow  resistance  across  the 
electrode.  The  electrolyte  has  several  order-of-magnitude  higher 
resistivity  than  the  electrodes  and  the  interconnect.  Thus,  in  an 
SOFC,  the  ohmic  polarization  is  typically  dominated  by  ion  resis¬ 
tance  thorough  the  electrolyte.  The  ohmic  overpotential  can  be 
expressed  as  [23]: 


2.2.2.  Cathode  activation  overpotential 

The  overall  oxygen  reduction  and  incorporation  at  the 
electrode-electrolyte  interface  can  be  written  as 

!o2(g)  +  Vo-  (el)  +  2 e-(c)  4*  0*(el) 

where  V0##(el)  and  0^(e/)  denote  the  oxygen  vacancies  and  lattice 
oxygen  ions  in  the  bulk  of  the  electrolyte  and  e-(c)  are  the  elec¬ 
trons  within  the  cathode.  As  with  the  oxidation  of  hydrogen  at  the 
anode,  the  global  reaction  is  the  result  of  elementary  steps.  Here,  it 
is  assumed  that  oxygen  reduction  proceeds  in  two  steps  [9]: 

•  Adsorption/desorption: 

02(g)  +  2(c)  ^  20ad(c)  (R6) 


hohm  ~  ie^el 

(31) 

Re,  =  =(i2cm2) 

^el 

(32) 

ael  =  croP"1  exp  (-^)  =  (Scm-1) 

(33) 

3.  Simulation  procedure 

The  model  described  in  Section  2  is  used  to  determine  the  cell 
voltage  for  a  given  current  density,  from  zero  current  to  the  limiting 
current.  The  cell  potential  is  expressed  as  the  difference  between 
the  equilibrium  potential  Erev  and  the  sum  of  all  the  relevant  over¬ 
potentials,  which  depend  on  the  current  density. 
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The  solution  proceeds  as  follows. 


(1 )  We  calculate  the  equilibrium  potential  based  on  the  global  elec¬ 
trochemical  reaction  at  zero  current. 

The  equilibrium  potential  depends  on  the  fuel  and  oxidant 
compositions  in  the  corresponding  channels,  temperature  and 
pressure. 

(2)  We  calculate  the  concentration  overpotentials. 

•  First  we  set  up  the  boundary  conditions: 

The  boundary  conditions  at  the  channel-electrode  inter¬ 
face  are  established  by  requiring  the  gas-phase  species 
concentrations  to  match  those  in  the  gas  channel.  At  the 
anode/electrolyte  interface,  the  boundary  condition  of  each 
gas  species  depends  on  the  charge-transfer  chemistry,  and 
hence  the  current  density. 


Jk  =  vk±  (34) 

•  We  calculate  the  molar  flux,  Jk  using  the  current  values  of 
the  concentrations  by  substituting  in  the  DGM,  which  is  rear¬ 
ranged  in  the  matrix  form  as  follows: 


-Vq  - 


Bo 


uv  De 

r  mix  uiM 


yh' VP  =  E  -  W  +A  =  [H]lJ] 

K  '  >)  ij  M 

(35) 


where  [/]  is  a  flux  column  vector  and  [H]  matrix  is  defined  as 


(36) 


Because  the  pressure  is  determined  using  the  ideal  gas  law,  the 
left  hand  side  of  Eq.  (35)  is  calculated  from  the  present  values  of 
the  concentrations  of  gas  species  and  the  concentration  boundary 
conditions  at  the  channel,  estimating  vq  using  forward  difference 
approximation.  The  matrix  [H]  in  the  right  hand  side  of  Eq.  (35)  is 
expressed  in  terms  of  the  diffusion  coefficients  and  the  mole  frac¬ 
tions  in  Eq.  (36)  and  calculated  in  the  same  way.  Hence,  the  molar 
flux  can  be  calculated  as 

L/]  =  [H]-1[LHS]  (37) 

We  then  substitute  Jk  in  the  following  conservation  equation  of  the 
gas  species: 

fir 

-gf  =  AsSsurLk  ~  v  Jk  (k  =  1 .....  Kg)  (38) 

(k  =  1, . . . ,  Ks)  (39) 

The  conservation  equations,  Eq.  (38)  and  (39),  become  a  set  of  ordi¬ 
nary  differential  equations  when  Eq.  (37)  is  cast  into  a  finite-volume 
form  using  the  flux  boundary  conditions  at  the  interface.  The  elec¬ 
trodes  are  approximated  as  continuous  media,  with  homogenized, 
volume  averaged  properties. 


•  We  calculate  sk  using  the  heterogeneous  thermochemistry  model 
based  on  the  current  values  of  the  concentrations  of  the  gas  and 
surface  species. 

The  ODEs  of  the  gas  species  and  surface  species  are  solved 
simultaneously  using  ‘odel5s’  function  in  MATLAB®  with  the 
Gear’s  method  option  on. 

•  We  calculate  the  concentrations  at  the  interface  between  the 
electrodes  and  the  electrolyte  when  the  solution  reaches  steady 
state. 

•  We  calculate  the  concentration  overpotentials  for  each  electrode, 

Picon, a  ^nd  TJcon,c • 


Fig.  1.  Terms  in  the  DGM  that  have  been  modified  and  the  impact  of  the  changes  in 
the  coefficient  on  the  contribution  of  the  corresponding  terms  to  the  concentration 
gradient. 


(3)  We  calculate  the  activation  overpotentials  for  each  electrode: 

The  activation  overpotentials,  r]a,a  and  rja,c ,  are  computed 
from  BV  type  i-V  relations  for  the  anode  and  cathode,  respec¬ 
tively.  rja,a  and  r\a,c  are  determined  from  the  nonlinear  solver 
function,  ‘fsolve’  function  in  MATLAB®. 

(4)  The  ohmic  overpotential  is  calculated  using  Eq.  (33). 

(5)  The  cell  voltage  is  calculated  using  Eq.  (2). 

4.  Simulation  results 

Among  the  numerous  ID  MEA  models  developed  to  describe 
the  transport-chemistry  interaction  within  the  different  elements, 
the  model  of  Zhu  et  al.  [9]  is  the  most  detailed  one.  Our  model, 
while  using  the  same  equations,  modifies  the  transport  and  acti¬ 
vation  overpotential  models  to  improve  their  accuracy,  as  shown 
below.  In  order  to  evaluate  the  effect  of  each  improvement,  we  com¬ 
pute  i-V  curve  using  a  single  rate-limiting  reaction  mode  across  the 
entire  current  density  range,  similar  to  that  used  in  Zhu  et  al.  [9]. 
Thus,  Reaction  (R3)  is  assumed  to  be  the  rate-limiting  reaction  over 
the  entire  range  of  current  density.  The  button-cell  experimental 
results  by  Jiang  and  Virkar  [7]  are  used  to  establish  the  empirical 
parameters,  i ^  and  i*  2 ,  in  the  electrochemistry  model  by  fitting  the 
MEA  experimental  performance  when  mixtures  of  hydrogen  and 
steam  are  fed  to  the  anode.  First,  we  demonstrate  the  limitation  of 
this  model,  even  after  modifying  the  transport  coefficients.  Next, 
we  use  two  different  reactions  as  rate-limiting  reactions,  Reaction 
(R3)  at  low  current  density  and  Reaction  (R1 )  at  high  current  den¬ 
sity,  and  show  that  this  model  predicts  the  limiting  current  density 
more  accurately. 

4.1.  Model  with  a  single  rate-limiting  reaction 

First,  we  investigated  the  impact  of  the  corrected  transport  coef¬ 
ficients,  and  compared  the  results  with  those  obtained  in  [9].  We 
corrected  the  effective  Knudsen  diffusion  coefficient  based  on  [3] 
and  adopted  the  permeability  reported  in  [6].  The  effective  Knudsen 
diffusion  coefficient  used  in  [9]  is  a  factor  of  2  larger  than  reported 
in  [3].  Furthermore,  the  permeability  used  in  [9]  is  smaller  than  the 
permeability  calculated  using  the  Carman-Kozeny  model  shown  in 
Eq.  (19)  by  a  factor  equal  to  0.4r. 

The  values  of  the  corresponding  terms  impacted  by  these  coef¬ 
ficients  are  two  to  three  times  larger  than  those  reported  in  [9], 
based  on  a  tortuosity  of  3.5.  Each  term  that  has  been  modified  in 
our  model  is  underlined  in  Fig.  1  and  the  approximate  influence  of 
each  term  is  shown  in  the  boxes  below  the  term. 

To  assess  the  impact  of  the  corrected  Knudsen  diffusion  coeffi¬ 
cient  and  the  permeability  in  the  DGM  model,  the  same  fuel  cell 
used  in  [9]  is  modeled  with  the  same  fitting  parameters  in  Table  2 
reported  in  that  paper. 
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Table  2 

Simulation  parameters. 


Parameters 

Value 

Units 

Parameters  for  an  SOFC  MEA  structure 

Anode 

Thickness  ( la ) 

1220 

|jim 

Porosity  ( e ) 

0.35 

Tortuosity  (r) 

3.50 

Average  pore  diameter  (do) 

1.00 

pan 

Average  particle  diameter  (dp) 

2.50 

|jim 

Specific  catalyst  area  (As) 

1080 

cm-1 

i*2  (Reaction  (R3)  is  rate-limiting) 

8.5 

A  cm-2  atm_3/4 

i^2  (Reaction  (Rl)  is  rate-limiting) 

15.963 

A  cm-2  atm-1 

Cathode 

Thickness  (Zc) 

30 

fxm 

Porosity  (s) 

0.35 

Tortuosity  (r) 

3.50 

Average  pore  diameter  (do) 

1.00 

fxm 

Average  particle  diameter  (dp) 

2.50 

fxm 

*o2 

2.8 

A  cm-2  atm-1/2 

Electrolyte 

Thickness  (Ze/) 

25 

fxm 

Activation  energy  of  O2-  (Ee;) 

8.0  x  104 

J  mol-1 

Pre-factor  of  O2-  (cr0) 

3.6  x  105 

Scirr1 

Operating  conditions 

Pressure 

1 

atm 

Temperature 

1073 

K 

Fig.  2  shows  substantial  difference  between  our  model  predic¬ 
tions  using  the  corrected  coefficients  and  those  reported  in  [9], 
especially  at  the  higher  current  density  when  the  concentration 
overpotential  is  most  relevant.  For  three  different  fuel  concentra¬ 
tions  at  the  anode,  our  model  reaches  the  limiting  current  density 
faster  than  that  reported  in  [  9  ] ,  and  the  difference  is  larger  for  higher 
fuel  concentrations. 

These  results  can  be  explained  as  follows.  As  the  current  density 
increases,  the  molar  flux  and  the  pressure  gradient  become  larger 
leading  to  higher  values  of  the  magnitudes  of  the  J  term  and  the 
B0  term  in  Fig.  1,  making  the  impact  of  the  two  coefficients  more 
important.  When  we  calculated  the  J  term  and  B  0  term  while  chang¬ 
ing  the  current  density  from  0.2  to  3  A  cm-2  for  the  case  where  the 
hydrogen  composition  is  50%,  the  values  of  these  two  terms  were 


Fig.  2.  Comparison  between  the  polarization  curves  computed  using  our  model  and 
those  reported  in  [9]. 


Fig.  3.  Comparison  between  the  experimental  data  [7]  results  shown  by  a  contin¬ 
uous  line,  our  model  results  shown  by  open  circles  and  Zhu  et  al.  [9]  model  results 
shown  by  a  broken  line  when  the  anode  stream  is  50%  H2  +  50%  H20. 

found  to  increase  by  an  order  of  magnitude.  Most  of  the  difference 
between  our  model  results  and  those  reported  in  [9]  comes  from 
the  J  term  which  includes  the  Knudsen  diffusion  coefficient.  On  the 
other  hand,  the  B0  term  is  not  as  significant  because  the  perme¬ 
ability  is  small,  O(10-15),  even  though  pressure  gradient  is  high, 
O(106). 

Figs.  3  and  4  show  that  correcting  the  transport  parameters  leads 
to  a  better  match  with  the  experimental  results  of  Jiang  and  Virkar 
[7].  Flowever,  even  with  these  corrections,  the  model  fails  to  accu¬ 
rately  predict  the  value  of  the  limiting  current  density.  In  the  next 
section,  we  show  that  using  a  different  reaction  step  at  high  current 
densities  can  further  improve  the  model  prediction. 


Fig.  4.  Comparison  between  our  model  results,  shown  by  lines,  and  experimental 
results  of  Jiang  and  Virkar  [7],  shown  by  open  circles,  rectangles,  and  triangles. 
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Fig.  5.  The  anode  activation  overpotential;  the  continuous  line  corresponds  to  the 
case  when  Reaction  (R1 )  is  rate-limiting  and  the  broken  line  corresponds  to  the  case 
when  Reaction  (R3)  is  rate-limiting  while  the  anode  stream  is  50%  H2  +  50%  H20. 


Fig.  6.  Comparison  between  our  model  and  experimental  results  by  Jiang  and  Virkar 
[7]. 


4.2.  Model  with  the  rate-limiting  reaction  switch-over 

Arguments,  shown  in  Section  2.2.1,  supporting  the  use  of  Reac¬ 
tion  (R3)  as  the  rate-limiting  reaction  may  not  be  valid  near  the 
limiting  current  density.  It  has  been  claimed,  based  on  exper¬ 
imental  results  in  [16,21],  that  the  adsorption  of  hydrogen  at 
the  anode  may  be  rate-limiting  process.  We  have  also  presented 
arguments  to  the  effect  that  hydrogen  adsorption  can  be  the  rate- 
limiting  reaction  under  the  condition  of  high  current  density.  At 
high  current  density,  a  vast  amount  of  water  is  produced  and  may 
compete  against  hydrogen  for  adsorption  sites,  reducing  hydrogen 
adsorption  rate  [1].  Therefore,  we  propose  a  new  model  in  which 
the  rate-limiting  reaction  switches  from  the  hydroxyl  oxidation 
charge-transfer  reaction,  Reaction  (R3),  to  the  hydrogen  adsorption 
reaction,  Reaction  (R1 ),  near  the  limiting  current  density. 

Fig.  5  shows  the  anode  activation  overpotential  when  Reaction 
(Rl)  or  (R3)  is  assumed  to  be  the  rate-limiting  reaction  over  the 
entire  range  of  current  density.  Results  show  that  at  lower  current 
densities,  Reaction  (R3)  overpotential  is  indeed  the  larger  value 
because  it  is  the  more  sluggish  reaction.  However,  as  the  current 
density  increases,  the  concentration  of  hydrogen  at  the  interface 
decreases  resulting  in  the  reduction  of  the  exchange  current  den¬ 
sity  and  Reaction  (Rl )  overpotential  rises  sharply.  This  is  consistent 
with  the  form  of  the  BV  equations  shown  in  Table  1 :  the  power  of 
the  hydrogen  concentration  in  i0  is  1  when  Reaction  (Rl )  is  assumed 
to  be  the  rate-limiting  reaction,  while  it  is  1  /4  when  Reaction  (R3)  is 
assumed  to  be  the  rate-limiting  reaction.  When  hydrogen  adsorp¬ 
tion  is  rate-limiting,  the  anode  activation  overpotential  soars  as  the 
exchange  current  density  decreases  because  of  the  drop  in  hydro¬ 
gen  concentration. 

In  the  following  calculation,  we  use  the  point  of  intersection 
of  the  two  curves  in  Fig.  5  as  the  switch-over  point  between  two 
reactions,  that  is,  the  point  at  which  we  assume  that  Reaction  (R3) 
is  no  longer  the  rate-limiting  reaction,  and  instead  Reaction  (Rl )  is 
the  rate-limiting  reaction. 

Results  of  our  extended  model  are  shown  in  Fig.  6.  The  results  of 
the  extended  model  agree  better  with  the  experimental  measure¬ 
ments,  although  the  limiting  current  density  is  now  slightly  under 
predicted,  especially  in  the  case  of  higher  hydrogen  concentration. 
This  can  be  explained  as  follows.  When  Reaction  (Rl )  was  assumed 
to  be  the  rate-limiting  reaction,  the  numerical  value  of  i*  2  in  Table  1 


was  obtained  by  fitting  the  experimental  results  for  the  case  of  34% 
H2  in  the  anode  stream  (and  66%  H20).  This  value  of  zj^  was  then 
assumed  to  be  independent  of  the  hydrogen  concentration  and  was 
used  in  all  the  other  cases,  i.e.  for  20%  and  50%  H2  in  the  anode 
stream.  Note  that  ifr  p h2  is  the  hydrogen  adsorption  rate.  Assuming 
that  i^2  is  constant  is  consistent  with  the  assumption  in  Section 
2.2.1  that  only  hydrogen  can  adsorb  on  the  Ni  surface.  However, 
when  we  consider  the  competitive  adsorption  on  the  Ni  surface,  the 
TPB  may  be  covered  with  more  water  molecules  when  the  water 
concentration  in  the  anode  stream  is  higher.  Thus,  zj^  may  depend 
on  the  concentration  of  H20  and  its  value  may  be  higher,  instead 
of  being  constant,  when  there  is  50%  H20  in  the  anode  stream  than 
when  there  is  66%  H20.  In  other  words,  the  activation  overpotential 
may  be  over  predicted  in  the  case  of  50%  H2  when  the  value  of  z'j^  is 
estimated  using  the  experimental  data  with  34%  H2.  This  explains 
why  the  limiting  current  density  is  slightly  under  predicted  in  Fig.  6 
when  the  anode  stream  is  50%  H2. 

When  the  hydrogen  adsorption  at  the  anode  is  taken  to  be 
the  rate-limiting  reaction,  the  exchange  current  density  corre¬ 
sponds  directly  to  the  current  density  of  the  cell  if  all  the  adsorbed 
hydrogen  is  consumed  in  the  electrochemical  reaction,  as  shown 
in  Table  1.  In  other  words,  the  exchange  current  density  is  the 
maximum  possible  current  density  when  the  hydrogen  adsorp¬ 
tion  reaction  is  rate-limiting.  This  is  consistent  with  the  i-r)act,a 
relationship  when  Reaction  (Rl)  is  assumed  to  be  rate-limiting: 


i  =  io 


1  -  exp 


^  2Fr)act,q^j 


(40) 


As  shown  in  Eq.  (40),  a  significantly  higher  activation  overpotential 
is  required  to  support  a  current  density  close  to  the  exchange  cur¬ 
rent  density.  Therefore,  the  limiting  current  density  is  very  close  to 
the  exchange  current  density. 


Table  3 

Rate-limiting  switch-over  point. 
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Fig.  7.  The  contributions  of  the  five  overpotentials  when  the  anode  stream  is  50% 
H2  +  50%  H20. 


Table  3  shows  the  computed  limiting  current  density,  and 
the  switching  current  density,  iSWitch>  that  is  the  current  density  at 
which  we  switch  between  Reactions  (Rl)  and  (R3).  Results  show 
that  iSWitch  =  Q-6i -o*  Calculations  in  Table  3  show  that  the  switching 
current  density  is  nearly  80%  of  the  limiting  current  density. 

We  conclude  this  section  by  restating  the  result  of  using  dif¬ 
ferent  reactions  as  rate-limiting  in  different  ranges  of  the  current 
density.  When  the  hydrogen  adsorption  rate  is  sufficiently  fast  to 
support  the  current  density,  such  as  1.67  times  the  current  den¬ 
sity,  the  anode  activation  overpotential  is  determined  by  Reaction 
(R3).  As  the  current  density  increases  and  the  corresponding  hydro¬ 
gen  adsorption  rate  decreases  below  this  level,  the  rate-limiting 
reaction  switches  to  Reaction  (Rl).  In  this  regime,  the  anode  acti¬ 
vation  overpotential  rises  abruptly  with  the  current  density  and  the 
electrochemical  reaction  proceeds  very  slowly  until  it  is  no  longer 
possible. 


0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8 

Current  density  (A  cm'2) 

Fig.  8.  The  contributions  of  the  five  overpotentials  when  the  anode  stream  is  34% 
H2  +  66%  H20. 


Current  density  (A  cm-2) 


Fig.  9.  The  contributions  of  the  five  overpotentials  when  the  anode  stream  is  20% 
H2  +  80%  H20. 

4.3.  Contribution  of  each  overpotential 

Figs.  7-9  show  the  contributions  of  each  overpotential  when 
the  hydrogen  concentration  at  the  anode  gas  channel  is  50%,  34%, 
and  20%,  respectively.  When  the  hydrogen  composition  is  50%  or 
34%,  the  ohmic  overpotential  and  the  cathode  activation  overpo¬ 
tential  have  almost  the  same  magnitude,  which  is  twice  as  large 
as  the  anode  concentration  overpotential.  They  are  also  more  than 
twice  as  large  as  the  anode  activation  overpotential,  except  near 
the  limiting  current  density  where  anode  activation  overpotential 
is  significant.  When  the  fuel  composition  is  20%  H2,  all  overpoten¬ 
tials,  except  the  cathode  concentration  overpotential,  have  nearly 
the  same  magnitude  in  the  region  away  from  the  limiting  current 
density.  For  all  cases,  as  expected,  the  cathode  concentration  over¬ 
potential  is  negligible. 


Fig.  10.  Comparison  between  the  experimental  results  by  Jiang  and  Virkar  [7]  and 
the  calculated  operating  cell  voltages  corresponding  to  the  single  rate-limiting  reac¬ 
tion  mode  and  the  rate-limiting  reaction  switch-over  mode,  respectively,  when  the 
anode  stream  is  34%  H2  +  64%  H20. 
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Table  4 

The  hydrogen  concentrations  at  the  channel,  and  the  limiting  current  densities  and  the  corresponding  hydrogen  concentrations  at  the  interface  for  the  single  rate-limiting 
reaction  mode  and  the  rate-limiting  reaction  switch-over  mode. 


Channel 

Single  rate-limiting 

Rate-limiting  switch-over 

Mole  fraction  (%) 

Concentration 
(mol  m-3) 

Limiting  current 
density  (A cm-2) 

Concentration  @  interface 
(mol  m-3) 

Limiting  current  density 

(A  cm-2) 

Concentration  @  interface 
(mol  ur3) 

20 

2.17 

1.55 

0.002 

1.05 

0.738 

34 

3.74 

2.56 

0.002 

1.75 

1.256 

50 

5.53 

3.65 

0.017 

2.03 

1.817 

For  all  known  hydrogen/ oxygen  fuel  cells,  regardless  of  elec¬ 
trolyte  type,  it  is  believed  that  electrochemical  reactions  at  the 
cathode  are  rate  limiting  and  the  activation  overpotential  is  almost 
entirely  due  to  cathode  [24].  This  is  because  hydrogen  electro¬ 
oxidation  is  extremely  rapid  on  a  wide  range  of  catalysts.  However, 
our  results  show  that  the  ohmic  overpotential  has  the  same 
magnitude  as  the  cathode  activation  overpotential.  The  ohmic  over¬ 
potential  is  dominated  by  the  ionic  resistance  in  the  electrolyte, 
which  is  a  thermally  activated  vacancy  hopping  mechanism.  The 
relatively  lower  operating  temperature  of  800  °C  might  explain 
these  results.  As  we  expect,  the  cathode  activation  overpotential 
is  larger  than  the  anode  activation  overpotential. 

However,  we  also  see  that  the  electrochemical  reaction  can¬ 
not  be  sustained  near  the  limiting  current  density  because  of  the 
anode  activation  overpotential.  It  is  often  assumed  that  the  con¬ 
centration  overpotential  captures  the  impact  of  transport  and  the 
activation  overpotential  accounts  for  electrochemistry.  Thus,  it  has 
been  argued  that  the  limiting  current  density  occurs  because  of 
transport  limitation,  and  the  limiting  current  density  has  been 
explained  using  the  rapid  increase  of  the  concentration  overpoten¬ 
tial.  However,  the  activation  overpotential  also  reflects  the  impact 
of  the  transport  since  the  exchange  current  density  depends  on  the 
concentrations.  Fig.  10  shows  a  comparison  between  the  experi¬ 
mental  results  by  Jiang  and  Virkar  [7]  and  the  calculated  operating 
cell  voltages  using  the  single  rate-limiting  reaction  mode  and  the 
rate-limiting  reaction  switch-over  mode,  with  their  correspond¬ 
ing  dominant  overpotential.  For  the  single  rate-limiting  mode,  the 
limiting  current  density  is  governed  by  the  rapid  increase  in  the 
concentration  overpotential  when  the  hydrogen  concentration  at 
the  interface  approaches  almost  zero,  as  shown  in  Table  4.  In 
the  rate-limiting  reaction  switch-over  mode,  the  sharp  increase 
in  the  activation  overpotential  occurs  before  the  rapid  increase 
in  the  concentration  overpotential.  The  hydrogen  concentrations 
at  the  interface  are  not  low  enough  to  force  the  rapid  increase 
in  the  concentration  overpotential,  as  shown  Table  4.  However, 
the  current  density  is  close  enough  to  the  exchange  current  den¬ 
sity  to  cause  a  sharp  increase  in  the  activation  overpotential. 
Hence,  the  limiting  current  density  is  governed  mainly  by  the 
activation  overpotential  rather  than  the  concentration  overpoten¬ 
tial. 

5.  Conclusions 

The  model  presented  in  this  paper,  which  is  similar  in  construc¬ 
tion  to  the  model  presented  in  [9],  is  based  on  the  DGM,  detailed 
heterogeneous  thermo-chemistry  models,  and  detailed  electrode 
kinetics  models.  We  corrected  the  effective  Knudsen  diffusion 
coefficient  and  the  permeability  values,  and  showed  improved  pre¬ 
dictive  accuracy.  We  additionally  analyzed  the  possibilities  that 
a  different  intermediate  step  in  the  hydrogen  electrochemical 
oxidation  model  is  rate-limiting  and  proposed  the  rate-limiting 
switch-over  model.  The  new  model  substantially  improves  the  pre¬ 
diction  of  the  limiting  current  density  and  shows  better  match  with 
experimental  results. 


Still,  there  are  many  improvements  that  can  further  contribute 
to  the  model  accuracy.  A  uniform  temperature  is  imposed  through¬ 
out.  However,  non-uniform  heat  release  due  to  thermo-chemistry 
and  various  overpotential  losses  may  introduce  temperature  non¬ 
uniformity.  Note  that  the  steam  reforming  is  highly  endothermic, 
and  the  ohmic  resistance  associated  with  ion  transport  through 
the  electrolyte  and  the  electrochemical  reactions  are  exother¬ 
mic. 

Recently,  it  was  reported  that  adding  an  interlayer  between  the 
electrode  and  the  electrolyte  improves  the  performance  of  SOFC 
[25].  Even  though  Zhao  and  Virkar  [26]  analyzed  the  impact  of  the 
interlayer,  more  work  is  needed  to  explain  the  role  of  the  interlayer. 

Although  our  anode  transport  and  thermo-chemistry  models 
assumed  that  the  fuel  is  methane  and  the  anode  is  Ni/YSZ,  it 
was  validated  only  for  hydrogen.  Using  methane  as  a  fuel,  it  has 
been  observed  that  SOFC  is  rapidly  deactivated  due  to  carbon 
deposition  on  the  anode  because  nickel  in  the  anode  catalyzes 
the  formation  of  graphite  from  hydrocarbons  and  its  deposi¬ 
tion  on  the  surface.  However,  almost  all  experimental  results 
of  SOFC  using  methane  as  a  fuel  were  conducted  while  using 
copper  in  the  anode  material  [27-30]  to  prevent  carbon  depo¬ 
sition.  Those  results  with  a  different  catalyst,  copper,  prevent 
us  from  validating  our  model  applying  the  detailed  heteroge¬ 
neous  kinetics  of  methane  on  the  nickel  surface.  The  detailed 
heterogeneous  kinetics  of  methane  on  the  copper  surface  may  be 
valuable. 

The  concentrations  imposed  as  boundary  conditions  at  channels 
might  be  different  from  those  in  the  inlet  flow  stream  considering 
the  flow  pattern  inside  the  button  cell  experiment.  The  inaccurate 
concentration  boundary  condition  at  channels  might  contribute  to 
discrepancies  between  the  model  simulation  results  and  experi¬ 
mental  results. 
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